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Abstract 

We simulate lattice QCD with an irrelevant chiral 4-fermion interaction which allows us to 
simulate at zero quark mass. This enables us to study the finite-temperature chiral-symmetry- 
restoring phase transition for 2 massless quark flavours, which is believed to be second order. In 
particular, it enables us to estimate the critical exponents which characterize the universality class 
of this transition. Our earlier simulations on A't = 4 and Nf = 6 lattices revealed that finite lattice- 
spacing artifacts on such coarse lattices affect the nature of the transition. We are now simulating 
on Nt = 8 lattices (12^ x 8, 16^ x 8 and 24^ x 8 lattices) where we expect to expose the continuum 
behaviour of this transition. 
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I. INTRODUCTION 



In the standard formulations of lattice QCD, the Dirac operator becomes singular in the 
chiral (zero quark mass) limit. This prevents one from simulating at zero quark mass, and 
makes simulations at small quark masses prohibitively expensive. In addition, since the 
chiral condensate vanishes on a configuration by configuration basis in this limit, even if one 
were able to perform zero-mass simulations, such simulations would tell us nothing about 
spontaneous chiral symmetry breaking. 

We have modified the standard staggered- quark action for lattice QCD by the addition of 
a chiral 4-fermion interaction of the form of Gross-Neveu chiral extension of the Nambu- Jona- 
Lasinio model . The idea of combining gauge theories with 4-fermion interactions 

has been explored by others We choose our 4-fermion operator to preserve the 

reduced symmetries of the staggered-quark action. 

This new action, by generating a chiral-symmetry preserving "dynamical quark mass", 
renders the Dirac operator non-singular, even for zero quark mass. In addition, the chiral 
condensate also remains finite in the chiral limit on each configuration. The chiral con- 
densate vanishes as it must on any finite lattice, because the orientation of the condensate 
rotates from configuration to configuration such that the average over the infinite ensemble 
of configurations vanishes. This method of enforcing the absence of spontaneous symmetry 
breaking on a finite lattice is identical to that seen in simple spin models, giving us con- 
fidence that we should be able to use this action to study critical behaviour of QCD, in 
the same way that these spin models allow determination of their critical behaviour using 
Monte-Carlo methods (see for example ^). 

Two-fiavour QCD is expected to have a second-order finite temperature chiral-symmetry- 
restoring (deconfining) phase transition at zero quark mass. This transition is expected to lie 
in the universality class of the 3-dimensional 0(4) sigma model. With the reduced symmetry 
of the staggered lattice we would expect this transition to lie in the universality class of the 
0(2) sigma model, until the lattice spacing is fine enough that the transition is sensitive to 
the continuum symmetry. At non-zero quark mass, this transition is expected to soften to 
a mere crossover. Hence to study this second order transition, it is necessary to simulate at 
zero quark mass (or so close to zero that the chiral extrapolation is unambiguous). 

Our earlier simulations on A^^ = 4 and Nt = 6 lattices P, 0] showed the nature 
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of the transition to be sensitive to lattice artifacts on such coarse lattices. The transition 
for Nt = 4 could be made first or second order, depending on the size of the 4-fermion 
coupling. The Nf = 6 transition appeared to be second order, but the measured critical 
exponents were those characteristic of a tricritical point, rather than those of the 0{2)/0{4) 
spin models. We are now extending these simulations to A^t = 8 where the lattice artifacts 
will be appreciably smaller. 

In section 2 we define our new, "xQCD", staggered action. Section 3 describes our 
ongoing Nt = 8 simulations and presents some preliminary results. In section 4 we give our 
conclusions and plans for the future. 



II. LATTICE QCD WITH CHIRAL 4-FERMION INTERACTIONS 

The standard staggered quark action differs from the continuum action by terms (9(a^), 
where a is the lattice spacing. This means we have the freedom to add terms which only 
contribute at 0{a'^) relative to the standard action, i.e. dimension 6 or higher operators. All 
such additions will do is change the size of the O(a^) corrections and hence the approach to 
the continuum. However, such actions will have the same continuum limit as the standard 
action. 

Sllfl, III to the 



We choose to add a chiral 4-fermion interaction of the Gross-Neveu form 
standard staggered action. Since we are interested in 2-flavour QCD, we choose an action 
whose naive continuum form would be 

Cf = i,{Ip+ m)i, - _ (^75x37^)^]. (1) 

where A'^2x2 is the number of SU{2) x SU{2) flavour doublets. Since the fermion doubling 
in the staggered fermion lattice implementation leads naturally to 4 flavours, the actual 
continuum action will be 

Cf = m)i, - - (^75^5^)2]. (2) 

DiV4x4 

where A^4x4 is the number of SU (4) x SU (4) flavour quartets. Clearly these two forms are 
equivalent if we replace A^/2 since A^4x4 = N2X2/2. 

Now we turn to the staggered lattice implementation of QCD with such 4-fermion inter- 
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actions The fermion part of the action is 

^/ = E + m + 1 Y.i^. + ten,)]x + E ^Nni'y' + vr^), (3) 

sites i s ° 

where e{x) = (— 1)^+?'+^+*^ cr and vr are auxihary fields residing on the dual lattice and 
i runs over the 16 sites on the dual lattice closest to the site x. Integrating out these 
auxiliary fields gives the lattice version of equation |2l remembering that the total number 
of flavours Nf = AN^y.^. Here 7 = Since the fermion determinant becomes real in the 

continuum limit, we use the determinant of M^M for our simulations and take the 4th root 
of this determinant using the hybrid molecular-dynamics method with 'noisy' fermions to 
get Nf = 2. We refer to lattice QCD with this action "xQCD". 

It is easy to see that this action preserves the same U{1) vestige of the S'[/(4) x SU{4:) 
continuum chiral symmetry that is preserved by the standard staggered action at m = 0, 
which is why this particular form was chosen. 

The reason it is possible to simulate this action at zero quark mass is clear in the confined 
phase. Here ipip develops an expectation value. Hence since 

(a) = {i,^)h (4) 

the quark field develops an effective mass m = (a), which makes the Dirac operator non- 
singular. What is less clear is why the Dirac operator remains non-singular in the chirally 
symmetric phase. 

7 should be chosen large so that the 4-fermion coupling is small and thus any additional 
flavour symmetry breaking over that already present in the standard staggered action, is 
also small. In particular 7 should be chosen far larger than that for which there is a bulk 
transition in the absence of the gauge flelds. In addition we shall try to keep 7 large enough 
that the critical value of /5 = 6/(7^ is not too far from that without the 4-fermion term, 
i.e. the gauge interactions should be primarily responsible for quark binding and chiral 
symmetry breaking. 

III. "xQCD" AT FINITE TEMPERATURE. 

Since Nt = 4 and Nt = 6 simulations showed evidence that the nature of the flnite 
temperature phase transition was being affected by flnite lattice spacing artifacts, we are 
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now simulating at Nt = 8. In particular we are running on 12^ x 8, 16^ x 8 and 24^ x 8 
lattices. 

We are currently running at zero quark mass with 7 = 10. At this 7 we are able to use a 
molecular-dynamics time increment of dt = 0.05, and find that at the worst it takes 600-700 
conjugate gradient iterations to invert the Dirac operator. However, because we now have 
a true critical point in the infinite volume limit, we have to face critical slowing down. For 
this reason we ran for 50,000 time units for each j3 in the transition region and are extending 
these runs to 100,000 time units. 

Let us now look at our preliminary results. In figure^we show the Wilson Line (Polyakov 
Loop) and the chiral condensate as functions of (3 across the transition region. The chiral 
condensate plotted here is really \J {tptp)'^ — ('^Ts^sV')^- What we see is that the transition 
region, 5.530 ^ (3 ^ 5.545, is very narrow, which explains why it is so difficult to make 
quantitative measurements for standard actions at non-zero quark masses. The crossover 
appears quite smooth which suggests that it is second order. Long-time correlations, which 
are apparent when one looks at the time histories of the measured observables close to the 
transition, also point to the transition being second order. 

The chiral order parameter ^ipipY or '(c^)' is non-zero on each configuration, but its ori- 
entation in the {{ipip) ,i{ip'^5^^ip)) or equivalently ((o"), (vr)) plane rotates from configuration 
to configuration (here the () represents the lattice average for a given configuration). It is 
this rotation that makes the ensemble average zero as it must be for a finite lattice. In 
figure 121 we show the distribution of values of the chiral condensate during runs at several 
/3s on a 24'^ x 8 lattice. What we notice for (3 = 5.5325 and (3 = 5.535 is that these points 
are concentrated in an annulus of non-zero radius. This is an indication that these f3s lie 
in the low-temperature phase where the chiral order parameter is non-zero. For jS = 5.5375 
there is no such annulus and the distribution peaks near the origin. This suggests that 
f5 = 5.5375 lies in the high temperature phase, where the chiral order parameter is zero and 
chiral symmetry is restored. From these observations, we would predict that the critical /3, 
Pc lies in the range 5.5350 ^ /3c ^ 5.5375. We also note that these distributions appear to 
have a circular symmetry in the plane. This ensures that the ensemble average will be zero. 
To see the distribution in the radial coordinate more clearly we have histogrammed these 
distributions to equal-size bins in (ipip)'^ — which yields a quantity proportional 

to the density in the planes of figure |21 These histograms are presented in figure El 
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FIG. 1: Wilson Line and chiral condensates as functions of /3. 

Now let us examine the chiral condensate in more detail. In figure ID we plot the chiral 
condensates from figure Q Because of the way in which the chiral condensate rotates from 
configuration to configuration, any methods used to derive the ensemble averaged condensate 
which has the correct continuum limit in the ordered phase will, of necessity, smooth out 
the transition. It is then a matter of 'experimentation' to determine which method is easiest 
to extrapolate to the infinite volume limit. The figure shows the average of the magnitude 
of the condensate, which has enjoyed some success with analyzing spin models 0,^^. The 
other method which we have used is to average the squares of the condensates. 

To determine the first of the critical exponents we need to fit the chiral condensate 
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24^x8 LATTICE |S=5.5325 7=10 m=0 
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FIG. 2: Distribution of chiral condensates for a) /3 = 5.5325, b) /3 = 5.535 and c) /? = 5.5375. 



to the scaling form 



(#) = c(/3 - /3, 



(5) 



for P close to Pc- Pm serves to distinguish the universality class of the transition. Pm = 
0.384(5) for 0(4), Pm = 0.35(llfor 0(Tl = 1/4 for a tricritical point and p^ = 1/2 



for a mean-field transition 
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16 



18l |. While it is doubtful if we will accumulate 



enough statistics to distinguish 0(2) from 0(4) behaviour, we should be able to distinguish 
0(2)/0(4) behaviour from mean- field or tricritical behaviour. 

Fitting our current 24'^ x 8 'data' to this scaling form over the range 5.5275 < P < 5.5375 
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24^x8 LATTICE /S=5.5325 7=10 m=0 
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FIG. 3: Histograms of distributions of chiral condensates for a) /? = 5.5325, b) /? = 5.535 and c) 
/3 = 5.5375. 

gives (3rn = 0.4(1). This should improve when we finish extending our statistics from 50,000 
to 100,000 trajectories for all jSs in the scaling window. In addition, we have yet to extract the 
maximum information from our 'data' by using Ferrenberg-Swendsen interpolation methods 
[igj l to increase the number of f3 values contributing to our fit. 

Another method to obtain information on the position and nature of the transition is to 
use 4th-order Binder cumulants for our order parameter, for several lattice sizes. The 4-th 
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FIG. 4: Chiral condensates as functions of /3. 



order cumulant for the auxiliary fields is defined as 

-D4 ^ 



(6) 



where a = (a, vr). Here we use the auxiliary fields rather than ((V^V'), ^(V'Ts^sV')) because 
we only have a single stochastic estimator for these fermionic order parameters for each 
trajectory, and we would need at least 4 such estimators/configuration to be able to obtain 
an unbiased estimator of -B4. Note that B4 uses quantities measured on finite lattices, 
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and does not suffer from the difficulties one encounters with trying to estimate the chiral 
condensate itself. 

^4 = 1.242(2) for Q. = 1.460998... for a tricritical point and = 7t/2 for 

a mean-field transition [21]. i?4 = 2 for a crossover. i?4 = 1 for a first-order transition. 
These are all infinite volume results. While the B4 values associated with a crossover or a 
first-order transition are approached slowly with increasing volume close to the second-order 
transition, the infinite volume value is approached rapidly with increasing volume, at the 
transition. As /? ^ 0, B4 approaches its first-order value, while as (3 ^ 00 B4 approaches its 
crossover value. The rate at which these values are approached as one goes away from the 
transition increases with increasing volume. The curves for different volumes will cross very 
close to the phase transition, yielding both estimates for i?4 and hence of the universality 
class of this critical point, and of jSc for this transition. Our preliminary results for all 3 
lattice sizes are given in figure Although these measurements are very noisy, the curves 
for the 3 lattice sizes appear to cross at 5.5325 ^ /9c ^ 5.5375. The value at which they 
cross is uncertain enough to accommodate both 0(2) and tricritical behaviour. 

Perhaps the best way to estimate the position of the phase transition is from the peaks 
of the susceptibilities. For any observable X the corresponding susceptibility is defined by 

X^ = V{{X-{X)r). (7) 

We use Ferrenberg-Swendsen interpolation to locate the peak of the susceptibilities for the 
plaquette, the Wilson line and a. The advantage of this for determining /5c is that these 
peaks give estimates for the position of the transition on a finite volume, where it is merely 
a crossover. These jSc estimates are relatively insensitive to the lattice volume, which will 
make any attempt to extrapolate to infinite volume somewhat easier. For the 24^ x 8 lattice 
we estimate Pc = 5.535(1) while for the 16'^ x 8 lattice we get Pc = 5.537(2), by this method. 

IV. CONCLUSIONS 

Lattice QCD, where the fermion action is modified to include a chiral 4-fermion term 
shows promise for determining the precise nature of the Nf = 2 finite temperature phase 
transition for massless quarks, by determining the critical exponents of this transition. The 
ability to work at and near m = seems essential for extracting this critical behaviour. With 
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FIG. 5: Binder cumulants as functions of f3. 



conventional actions, the penalties for running at quark masses small enough to predict the 
chiral limit appear insurmountable, at least in the near future. To see this one needs only 
to look at the best 'data' currently available from simulations using state-of-the-art actions 



which do not permit using m = [2^ . We have already seen that simulations with A^^ = 4, 6 
lattices are plagued by finite lattice artifacts as was implicit in the analyses using standard 
actions. We are therefore using lattices with Nt = 8. 

Preliminary fits to the 24^ x 8 data yield critical exponent /3m = 0.4(1), consistent with 
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the expected 0{2) or 0(4) values. However, with an error this large, one cannot preclude 
tricritical or mean-field behaviour. More statistics are needed as well as reliable extrapola- 
tions to infinite volume. Use of Ferrenberg-Swendsen interpolation, thus making more use 
of the limited 'data' wc do have, is one way of effectively increasing our statistics. 

The scaling window is very narrow Af3 ^ 0.01, which is the principal reason why extrapo- 
lating from finite quark mass is extremely difficult. We are able to pinpoint the critical point 
from the peaks in the susceptibihties for the various observables. This method yields the 
best estimate of this position for a given lattice size and does not rely on some potentially 
unreliable infinite volume extrapolation. Since the volume dependence of /3c determined this 
way is minimal, this is probably also the best method for determining (3(. in the infinite 
volume limit. 

Binder cumulants show promise for determining the position and nature of the phase 
transition. However, these require much more statistics than we have obtained and possibly 
more than we can hope to obtain. They have the advantage that they make use of the chiral 
condensate as measured on the finite lattice, despite the fact that it averages to zero at all 
/?s. Such are the advantages of fiuctuation quantities. 

We will ultimately measure other critical exponents. One such, namely 6, measures the 
mass dependence of the condensate at /3c. Hence this will require simulations at small but 
finite quark mass. 

We are currently performing simulations with this action on zero temperature lattices. 
These simulations will calculate the spectrum of light hadrons, which should indicate the 
utility of this action for problems other than QCD thermodynamics, for which the ability 
to simulate at zero or physical quark masses is also important. 
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